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Abstract 

A continuum damage model for the prediction of damage onset and structural col- 
lapse of structures manufactured in fiber-reinforced plastic laminates is proposed. 
The principal damage mechanisms occurring in the longitudinal and transverse di- 
rections of a ply are represented by a damage tensor that is fixed in space. Crack 
closure under load reversal effects are taken into account using damage variables 
established as a function of the sign of the components of the stress tensor. Damage 
activation functions based on the LaRC04 failure criteria are used to predict the 
different damage mechanisms occurring at the ply level. The constitutive damage 
model is implemented in a finite element code. The objectivity of the numerical 
model is assured by regularizing the dissipated energy at a material point using 
Bazant’s Crack Band Model. To verify the accuracy of the approach, analyses of 
coupon specimens were performed, and the numerical predictions were compared 
with experimental data. 


Key words: Fracture Mechanics, Continuum Damage Mechanics, Composite 
Materials. 


1 INTRODUCTION 


The methodology for designing high-performance structures of composite ma- 
terials is still evolving. The complexity of the response of composite materials 
and the difficulty in predicting structural modes of failure result in the need 
for a well-planned test program. The recommended practice to mitigate the 
technological risks associated with such materials is to substantiate the per- 
formance and durability of the design in a sequence of steps known as the 



Building Block Approach (BBA)[1], The BBA ensures that cost and perfor- 
mance objectives are met by testing greater numbers of smaller less expensive 
specimens, assessing technology risks early in the program, and building on the 
knowledge acquired at a given level of structural complexity before progressing 
to a level of more complexity. 

Achieving substantiation of structural performance by testing alone can be 
prohibitively expensive because of the number of specimens and components 
required to characterize all loading conditions. BBA programs can achieve 
significant cost reductions by seeking a synergy between testing and analysis. 
The more the development relies on analysis, the less expensive it becomes. 

The use of advanced analytical or numerical models for the prediction of the 
mechanical behavior of composite structures can replace some of the mechan- 
ical tests and can significantly reduce the cost of designing with composites 
while providing to the engineers the information necessary to achieve an op- 
timized design. 

Strength-based failure criteria are commonly used to predict failure in compos- 
ite materials. A large number of continuum-based criteria have been derived to 
relate stresses and experimental measures of material strength to the onset of 
failure [2] -[4]. Failure criteria predict the onset of the different damage mech- 
anisms occurring in composites and, depending on the material, the geometry 
and the loading conditions, may also predict final structural collapse. 

For composite structures that can accumulate damage before structural col- 
lapse, the use of failure criteria is not sufficient to predict ultimate failure. 
Simplified models, such as the ply discount method, can be used to pre- 
dict ultimate failure, but they cannot represent with satisfactory accuracy 
the quasi-brittlc failure of laminates that results from the accumulation of 
several damage mechanisms. 

The study of the non-linear response of quasi-brittle materials due to the 
accumulation of damage is important because the rate and direction of dam- 
age propagation defines the damage tolerance of a structure and its eventual 
collapse. To model the phenomena of damage propagation, non-linear consti- 
tutive models defined in the context of the mechanics of continuum mediums 
have been developed and implemented in finite elements codes in recent years. 
The formalism of the thermodynamics of irreversible processes is a rigorous 
framework from which the constitutive models can be developed. 

The simplest way to describe damage is using a single scalar damage variable 
as proposed by Kachanov [5]. Damage can be interpreted as the creation of 
microcavities, and the damage variables as a measure of the effective surface 
density of the microdefects. Such a mechanical interpretation of damage as- 
sumes that the loads are resisted only by the undamaged ligaments in the 
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material. The stresses (a) in the ligaments, referred to as effective stresses, 
continue to increase until all ligaments are severed and the material has failed. 

The tensorial representation of damage is a formal and general procedure to 
represent the directionality of micro-cracks, which can take any direction in a 
medium depending on the load history, geometry, boundary conditions, and 
material properties. After Kachanov’s pioneering work, several damage models 
have been developed that describe damage as a second order tensor [6] -[9] or 
as a fourth order tensor [10]- [12], Second order tensors describe an initially 
isotropic material as an orthotropic one when damage evolves, whereas fourth 
order tensor models can remove all material symmetries and provide a more 
general procedure to simulate damage [13]. 

The application of continuum damage models in orthotropic or transversely 
isotropic materials, such as fiber-reinforced plastics (FRP), results in addi- 
tional difficulties. The nature and morphology of a material induces some pre- 
ferred directions for crack growth, i.e., crack orientations are not only induced 
by the loads, geometry and boundary conditions, but also by the morphology 
of the material. The interface between fiber and matrix is weaker than the 
surrounding material and interfacial debonding is normally the first damage 
mechanism to occur. Furthermore, residual thermal stresses occur in the com- 
posite plies due to different coefficients of thermal expansion of the fiber and 
matrix (micromechanical residual thermal stresses) and due to the different 
coefficients of thermal expansion in the longitudinal (fiber) and transverse 
(matrix) directions (macromechanical residual thermal stresses). 

Multiscale models and mesomodeling are two approaches used to evaluate 
the clastic and inelastic response of a material. Using homogenization laws, 
multiscale models define relations between a mesoscale, normally the scale 
of the finite elements, where material is considered homogeneous, and the 
microscale the scale of constituents, fiber and matrix. The constitutive models 
are defined at the microscale, and the strain and stress fields at the microscale 
and the mesoscale are related via transformation held tensors [14]- [18], or 
solved using finite elements [18], [19]. To reduce the amount of computations 
that need to be performed, periodicity of the material is invoked. 
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Mesomodeling is an alternative way to define damage models for composite 
materials that is more appropriate for large scale computations. Mesomodels 
treat the composite lamina [20]- [24] or sub-laminate [25] as a homogeneous 
material. When diffuse damage localizes in a narrow band and becomes a 
macro-crack, the response is dominated by the crack tip and its ability to dis- 
sipate energy. On the other hand, the material morphology, which is the main 
basis of homogenization techniques, loses importance due to the loss of peri- 
odicity. Therefore, in structures exhibiting stable crack propagation, i.e., when 
the macrocrack length does not increase under constant load, mesomodeling 
is more appropriate to predict the structural collapse than multiscale models. 

The main objective of the present paper is to develop a continuum dam- 
age model able to represent the quasi-brittle fracture of laminated composite 
structures, from damage onset up to final structural collapse. 

The majority of the material properties required by the present model can 
be measured using standard test methods. Most of the material properties 
that are required can be obtained from ply-based test methods. The use of 
ply properties, rather than laminate properties, is an advantage because it 
avoids the need to test laminates every time the lay-up or stacking sequence 
is modified. 

The proposed constitutive model accounts for crack closure under load reversal 
effects, an important phenomenon in cases where a composite structure is 
subjected to multiaxial loading. 

One important issue regarding the numerical modeling of damage is that the 
convergence of the solution through successive mesh refinement must be en- 
sured. The objectivity of the numerical model is ensured by adjusting the 
energy dissipated by each damage mechanism using a characteristic element 
length. The constitutive model proposed herein can be integrated explicitly, 
making it computationally efficient and, therefore, suitable to be used in large 
scale computations. 

This paper is organized as follows: a brief description of the damage mecha- 
nisms occurring in laminated composites is presented. Based on the mecha- 
nisms of damage identified, a new constitutive damage model is proposed. The 
constitutive model relates the damage mechanisms with a set of internal vari- 
ables. The constitutive model is implemented in a non-linear finite element 
code, and it is adjusted using a procedure based on Bazant’s Crack Band 
Model [26] that ensures the correct computation of the energy dissipated by 
each damage mechanism. The accuracy of the model is assessed by comparing 
the predictions with experimental data for an open-hole tension carbon-fiber 
specimen. 
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2 MECHANISMS OF DAMAGE AND FRACTURE IN LAMI- 
NATED COMPOSITES 


2. 1 Longitudinal failure 


In fiber-reinforced composites, the largest portion of the loads is resisted by the 
fibers. When these fail under either tension or compression, the internal loads 
must redistribute to other areas of the structure, and may cause a structural 
collapse. 

In composites with high fiber volume fraction and those where the strain to 
failure of the resin matrix is higher than the one of the reinforcing fiber, such as 
carbon-epoxy composites, longitudinal failures start by isolated fibre fractures 
in weak zones. The localized fractures increase the normal and interfacial shear 
stresses in the adjoining fibers, and the local stress concentrations promote 
matrix cracking, fibre matrix debonding and, for ductile matrices, conical shear 
failures [27]. When increasing the load further, additional fibre fractures occur, 
leading to final collapse. 

Failure under longitudinal tensile loading occurs in both constituents, and 
fracture occurs along a plane that is parallel to the fibers and the thickness 
direction (Figure 1). A simple non-interacting failure criterion based on max- 
imum stress or maximum strain along the longitudinal direction can usually 
provide an accurate measure of longitudinal tensile failure [3] . 



Fig. 1. Longitudinal failure in 0° plies [28]. 


Compressive failure of aligned fiber composites occurs from the collapse of 
the fibers as a result of shear kinking and damage of the supporting matrix 
[29], [30]. A kink band in a carbon-epoxy laminate resulting from compressive 
longitudinal stresses is shown in Figure 2. 
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0° ply 



0.1mm 

Fig. 2. Kink band in a 0° ply [31]. 

Argon [32] was the first to analyze the kinking phenomenon, based on the 
assumption of a local initial fiber misalignment. Fiber misalignment causes 
shear stresses between fibers that rotate the fibers, which in turn increase the 
shear stress and leads to instability. 

Recently, the calculation of the critical kinking stress has been significantly 
improved with a more complete understanding of the geometry of the kink 
band as well as the incorporation of friction and material nonlinearity in the 
analysis models [28], [33]. 


2.2 Transverse failure 


Failure in the transverse direction encompasses both matrix cracking and fiber- 
matrix debonding. Under the presence of transverse tensile stresses and in- 
plane shear stresses, the combined effect of small defects present in a ply 
such as small fiber-resin debonds, resin-rich regions, and resin voids, trigger a 
transverse crack that extends through the thickness of the ply (Figure 3). The 
transverse cracks are formed without disturbing the fibers: they occur at the 
fiber-resin interface and in the resin. 


90° ply 



Fig. 3. Transverse matrix crack in a 90° ply. 

A fundamental issue that needs to be considered is the effect of ply thickness 
on the ply strength, usually called the ’in-situ effect’. As shown in Parvizi [34] 
and Chang’s [35] experiments, the constraints imposed by the neighboring 
plies of different fiber orientations cause an apparent increase in the tensile 
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and shear strengths of a ply compared to those of an unconstrained ply. The 
scale effect for the in-plane shear strength is represented in Figure 4 , where it 
is apparent that the strength of a ply is a function of the number of 90° plies 
stacked together using the model proposed by Camanho et al. [36]. 

The in-situ effect is a deterministic size effect that can be represented using 
fracture mechanics models of plies containing defects [3], [37]. 



Fig. 4. In-situ shear strength of T300/1034-C CFRP. 

Experimental results have shown that moderate values of transverse compres- 
sion have a beneficial effect on the strength of a ply. This effect can be observed 
in the experimental results obtained by Swanson [38], [39] shown in Figure 5. 
The failure envelope calculated using the LaRC04 failure criterion [4] is also 
shown in Figure 5. 

When the in-plane shear stress is large compared to the transverse compressive 
stress, the fracture plane is perpendicular to the midplane of the ply. This 
mode of failure is referred to as Mode A in Figure 5. However, increasing 
the compressive transverse stress causes a change in the angle of the fracture 
plane, as shown in Figure 6. This mode of failure is referred to as Mode B 
in Figure 5. Normally, for carbon-epoxy and glass-epoxy composites loaded in 
pure transverse compression, the fracture plane is at an angle (fracture angle, 
ao) °f 53° ± 3° with respect to the thickness direction [40]. Therefore, matrix 
cracking does not occur in the plane of the maximum transverse shear stress 

(45”). 
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Fig. 5. Strength as a function of transverse compression and in-plane shear. 



0.1 mm 


Fig. 6. Matrix crack in a 45° ply created by in-plane compressive transverse stress 

[31]- 


2.3 Delamination 


Delamination is a common damage mechanism in multidirectional laminated 
composites due to their weakness in the thickness direction. The three di- 
mensional stress states that occur near geometric discontinuities such as ply 
drop-offs, stiffener terminations, flanges, bonded and bolted joints, and access 
holes promote delamination initiation. Delamination causes a reduction of the 
bending stiffness of a composite structure and, when compressive loads are 
present, promotes local buckling. 

Delamination models [41]- [44] always represent a discrete crack separating ele- 
ments as do, for concrete, the pioneering work of Dudgale-Barenblatt [45], [46]. 
The purpose of present work is to propose a Continuum Damage Mechanics 
(CDM) model for the calculation of the initiation and propagation of intralam- 
inar damage. Consequently, delamination damage is not considered here. 


3 CONSTITUTIVE DAMAGE MODEL 


The thermodynamics of irreversible processes is a general framework that can 
be used to formulate constitutive equations. It is a logical framework for incor- 
porating observations and experimental results and a set of rules for avoiding 
incompatibilities. In this section, we present a constitutive damage model for 
laminated composites that has its foundation in irreversible thermodynamics, 
and that uses the LaRC04 criteria as damage activation functions. 


3. 1 Complementary free energy and damage operator 


To establish a constitutive law, it is necessary to define a scalar function 
corresponding to the complementary free energy density in the material. This 
function must be positive definite, and it must be zero at the origin with 
respect to the free variables (the stresses) [47] . The proposed definition for the 
complementary free energy density is: 

r-, _ a \l , a 22 — rr rr -I ^12 i 

T ~ 2 (l — di) Ei ' 2 {1 — d 2 ) E 2 Ei 11 22 + 2 (1 - d 6 ) G 12 

+ (aqi< 7 ii + 0:22022) AT + (/?ncr 11 + /?22C r 22) AM ( 1 ) 

where £j, E 2 , zq 2 and G 12 are the in-plane clastic orthotropic properties of 
a unidirectional lamina. The subscript 1 denotes the longitudinal (fiber) di- 
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rection, and 2 denotes the transverse (matrix) direction. The damage variable 
d\ is associated with longitudinal (fiber) failure, whereas c /2 is the damage 
variable associated with transverse matrix cracking and d 6 is a damage vari- 
able influenced by longitudinal and transverse cracks, an and «22 are the 
coefficients of thermal expansion in the longitudinal and transverse directions, 
respectively. (3\ \ and P 22 are the coefficients of hygroscopic expansion in the 
longitudinal and transverse directions, respectively. AT and AM are the differ- 
ences of temperature and moisture content with respect to the corresponding 
reference values. 

The stress tensor a corresponds to the average stress tensor over a repre- 
sentative volume that is assumed to be much larger than the diameter of a 
fiber. 

To ensure the thermodynamically irreversibility of the damage process, the 
rate of change of the complementary free energy G minus the externally sup- 
plied work to the solid a : £ at constant strains, must not be negative: 


G — a : e > 0 (2) 

This inequality corresponds to the positiveness of the dissipated energy and 
has to be fulfilled by any constitutive model [47] . Expanding the inequality in 
terms of the stress tensor and damage variables gives: 


dG ■ 


( 3 ) 


Since the stresses are variables that can vary freely, the expression in the 
parenthesis must be equal to zero to ensure positive dissipation of mechanical 
energy. Therefore, the strain tensor is equal to the derivative of the comple- 
mentary free energy density with respect to the stress tensor: 


dG 



H : a + a AT + (3 AM 


( 4 ) 


The lamina compliance tensor can be represented in Voigt notation as: 


H 


d 2 G 
da 2 


1 ^21 

(1 — di) E\ E 2 

1 

E\ (1 — c/ 2 ) E 2 
0 0 


0 

0 

1 


( 5 ) 
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(1 — c/g) G\2 J 



The closure of transverse cracks under load reversal, also known as the unilat- 
eral effect, is taken into account by defining fonr damage variables associated 
with longitudinal and transverse damage. To determine the active damage 
variables, it is necessary to define the longitudinal and transverse damage 
modes as follows: 


d\ = 
C?2 = 


r L fgnf 

dl+ kill 
j (0-22) 
d2+ \<T22\ 


+ d\- 


(-^n) 

knl 


+ C?2- 


(—0-22) 

k22 1 


( 6 ) 


where ( x ) is the McCauley operator defined as (x) := (x + |x|) /2. 

The present model tracks damage caused by tension loads (d + ) separately 
from damage caused by compression loads (gL). Depending on the sign of the 
corresponding normal stress, a damage mode can be either active or passive. 

The model also assumes that the shear damage variable, d 6 , is not affected 
by the closure effect. Shear damage is caused mainly by transverse cracks and 
these do not close under shear stresses (cr 12 ). Transverse cracks are influenced 
by transverse stresses ( 022 ) producing the closure of cracks and a friction re- 
tention whereas longitudinal cracks produce the same effect under longitudinal 
stresses (<Tn) [48]. The effect of friction is neglected in the present model. 


3.2 Damage activation functions 


The determination of the domain of elastic response under complex stress 
states is an essential component of an accurate damage model. Based on the 
previously described mechanisms of crack generation in advanced composites, 
a strain space is considered where the material is linear elastic. In the present 
model, it is assumed that the clastic domain is enclosed by four surfaces, each 
of them accounting for one damage mechanism: longitudinal and transverse 
fracture under tension and compression. Those surfaces are formulated by the 
damage activation functions based on the LaRC03 and LaRC04 failure crite- 
ria. The LaRC03-04 failure criteria have been shown to represent accurately 
the physical process of damage onset in laminated composites. The LaRC04 
criteria represent an evolution of the LaRC03 criteria: some criteria such as 
the one for fiber kinking, are more accurate in LaRC04. However, the improve- 
ment in accuracy is associated with a significant increase in the computational 
effort. The present damage model uses a combination of both sets of criteria 
to achieve a compromise between accuracy and computational efficiency. The 
full details of the derivation and validation of the LaRC04 failure criteria are 
presented in references [3], [4]. 
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The four damage activation functions, ipv, associated with damage in the lon- 
gitudinal (TV = 1+,1— ) and transverse (TV = 2+,2— ) directions, are defined 
as: 


F 1+ = 0 i + - n+ < 0 ; Fi_ = (j)\— - n_ < 0 


( 7 ) 


F 2+ = 02+ - r 2+ < 0 ; F 2 - = 02- - r 2 - < 0 


where the loading functions 0jv (TV = 1+, 1— , 2+, 2— ) depend on the strain 
tensor and material constants (elastic and strength properties). The elastic 
domain thresholds (TV = 1+,1 — ,2+,2— ) take an initial value of 1 when 
the material is undamaged, and they increase with damage. The clastic domain 
thresholds are internal variables of the constitutive model, and are related to 
the damage variables cIm (TIT = 1+, 1— ,2+,2— , 6) by the damage evolution 
laws. The clastic domain threshold defines the level of elastic strains that can 
be attained before the accumulation of additional damage. 


3.2.1 Longitudinal tensile fracture 

The LaRC04 criterion for fiber tension failure is a non-interacting maximum 
allowable strain criterion defined as: 


Ei dn — ^12d22 

<Pi+ = = y 

yyj 1 y\_T 


( 8 ) 


where the effective stress tensor d is computed as a — H 0 1 : e. H 0 is the 
undamaged compliance tensor obtained from equation (5) using d\ = d 2 = 
d,Q = 0 . 


3.2.2 Longitudinal compressive fracture 

The LaRC03 failure criterion for longitudinal compressive fracture postulates 
that a kink band is triggered by the onset of damage in the supporting matrix. 
Under this circumstance, the fibres lose lateral support and fail under the 
effect of longitudinal compressive stresses. The initial fiber misalignment and 
the rotation of the fibers as a function of the applied stress state are the 
parameters used in the damage activation function. 

The damage activation function used to predict damage under longitudinal 
compression (dn < 0) and in-plane shear (fiber kinking) is established as a 
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function of the components of the stress tensor cr (m ) in a coordinate system 
(m) representing the fiber misalignment: 


L z-m 
22 


ra + v <?■ 

Sr 


(9) 


where the longitudinal friction coefficient can be approximated as [3]: 


L Sl cos (2a 0 ) 

V ~ 77 2 

r c cos^ a 0 

The components of the effective stress tensor in the coordinate system asso- 
ciated with the rotation of the fibers are calculated as: 



a 22 = &11 SIH 2 + °22 COS 2 pf — 2 |cr 12 | sin </? c cos ip c ^ 

d” 2 = (&22 — dn) sin ip c cos (p c + \d\ 2 \ (cos 2 (p c — sin 2 cp 0 ') 

where the absolute value of the shear stress is taken because the misalignment 
angle can be positive or negative. 

The misalignment angle (<p c ) is determined using standard shear and longi- 
tudinal compression strengths, Sl and Xq [3]: 



It should be noted that the LaRC03 failure criterion, derived to predict the 
onset of damage in laminated composites, calculates the misalignment angle 
as a function of the applied stress. The LaRC03 failure criterion is modified 
here assuming a constant misalignment angle, corresponding to the rotation of 
the fibres at failure under pure longitudinal compression. This modification of 
LaRC03 failure criterion assures that 4>i- is a monotonic increasing function 
under any state of proportional loading. 

It should be pointed out that two criteria are used in LaRC03 for fiber kinking: 
Equation (9) for < 7 ^, < 0 and a second equation for < 7 ^ > 0. However, the 
omission of the second equation results in a minor loss of accuracy because the 
equation is the same as equation (13. a) (below) with the stresses transformed 
into the misaligned coordinate frame. 
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3.2.3 Transverse fracture perpendicular to the laminate mid-plane (a o — 0 o j 

Transverse matrix cracks perpendicular to the mid-plane of the ply, i.e. with 
a o = 0° (Figures 6 and 5), are created by a combination of in-plane shear 
stresses and transverse tensile stresses, or in-plane shear stresses and small 
transverse compressive stresses. These conditions are represented by the fol- 
lowing LaRC04 failure criterion: 


02 + ~ 


(|+ 2 | + y+ 2 ) if O 22 < 0 


(13) 


where g is the fracture toughness ratio defined as: g = Glc 


Gnc ' 


3.2.4 Transverse fracture with cko = 53° 

The LaRC04 matrix failure criterion for transverse compressive stresses con- 
sists of a quadratic interaction between the effective shear stresses acting on 
the fracture plane: 


0 2 - 




if O ’ 22 < 0 


where the effective stresses r^ s and r e ^ are computed as [4]: 


(14) 


T e g = (-cr 2 2 COS (ck 0 ) (sin (a 0 ) - g T cos (a 0 ) cos ( 0 )) ^ 
rfs = (cos (a 0 ) ( 1 6-12 1 + V L °22 cos (a 0 ) sin ( 0 )) ^ 

The sliding angle 9 is calculated as [4]: 


(15) 


9 = arctan ( ^ ~ 1 2 ^ ) (16) 

V ^22 sin(Q'o) J 

The transverse shear strength and transverse friction coefficient can be ap- 
proximated as: 



Y c cos (ckq) sin (ct 0 ) + 


-1 

tan(2ao) 


cqs(q;o) 

tan(2ao) 


(17) 
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The fracture angle «o is approximately 53° in uniaxial compression. With 
increasing amounts of in-plane shear, the fracture angle diminishes up to about 
40° (Mode B in figure 5) and then abruptly switches to 0° (Mode A in figure 
5). To find the correct angle of fracture, a maximization of the LaRC03-04 
failure criteria as a function of a should be performed. However, in order to 
improve the computational efficiency of the present model, it is assumed that 
the fracture angle can only take one of two discrete values: 0° or 53°. 

The elastic domain in the <Tn, 022 , <J\i space represented by the LaRC04 failure 
criteria is shown in Figure 7. 



Fig. 7. Elastic domain in the an, < 722 , 012 space. 
3. 3 Dissipation 


The rate of energy dissipation per unit volume resulting from the evolution of 
damage is given by: 


“ — yri + wr^2 + — Yid\ + HdA + Yq(1q > 0 (18) 

ud 1 od 2 (yds 

The form of the complementary free energy defined in equation (1) assures 
that the thermodynamic forces (Ym) conjugated to their respective damage 
variables (o?m) are always positive: 
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(19) 


>4 = 


dG 

ddi 

dG 

dd,2 

dG 

dd 6 


2 

a n -> n 
2(1 -d 1 fE 1 ~ U 
2 

_22 N n 
2(l-d 2 ) 2 £ 2 - U 
2 

a ii > D 
2(l-d 6 ) 2 Gi2 - 


Therefore, the condition of positive evolution of damage variables {(1m > 0) is 
a sufficient condition for the fulfillment of the second law of thermodynamics. 

It is important to note that the proposed model does not generate spurious 
energy dissipation, i.e., the loss or gain of mechanical energy, under crack 
closure or opening. At load reversal, the time derivative of the damage variable 
is non-zero (< 1m ^ 0). Considering equation (18), the thermodynamical forces, 
Y m , associated with the damage variables, (1m , must be zero to avoid spurious 
energy dissipation at load reversal [49]. This condition is trivially satisfied in 
the present model (equation (19)). 

Damage evolution without energy dissipation is physically inadmissible. There- 
fore, it is necessary to avoid damage evolution when the corresponding con- 
jugated thermodynamic force is zero. Consider the load history represented 
in Figure 8: the material is loaded in transverse tension and shear to t\ and 
then loaded to £ 2 - At time £ 2 , the damage variable d 2 evolves because the corre- 
sponding damage activation function is activated. However, the corresponding 
thermodynamic force is zero {cr 22 = 0, Y 2 = 0). 



This non-physical response is avoided by modifying the longitudinal and trans- 
verse damage activation functions. The transverse damage activation function 
is modified using the following equation: 
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= nun 


( 20 ) 



022 I 

IT I 


where the constant £7 is equal to a 2 2 when 0 2 -=02+ (Figure 8). 

The longitudinal damage activation function is modified by taking into account 
that under shear dominated loads, matrix cracking is the first form of damage 
to occur. After matrix cracking, the transverse and shear stresses are zero and 
the fiber misalignment angle tp tends to 7t/4. Under these circumstances, the 
kink band criteria, equation (9), reads: 


<j> i_ = min 



+ ^22) 


Sl 



( 21 ) 


3-4 Damage evolution 


The evolution of the threshold values rjv is mathematically expressed by the 
Kuhn-Tucker conditions'. 


t'n A d-; igy < 0 ; r nFn — 0 (22) 

Neglecting viscous effects, the damage activation functions, equations (7), al- 
ways have to be non-positive. While the damage activation function F N is 
negative, the material response is elastic. When the strain state activates a 
criterion, Fjy = 0, it is necessary to evaluate the gradient 4>n- If the gradient 
is not positive, the state is one of unloading or neutral loading. If the gradient 
4>n is positive, there is damage evolution, and the consistency condition has 
to be satisfied: 


F n — O.x — r.v — 0 (23) 

Two important characteristics of the model proposed here are that the thresh- 
old values are a function of the damage variables, and that the loading func- 
tions depend on the strain tensor. Under these conditions, it is possible to 
explicitly integrate the constitutive model [11], [12], 

In the definition of the constitutive model, it is necessary to represent the 
relation between active and inactive clastic domains. The evolution of an active 
clastic domain is defined by the consistency condition, i.e., it is defined with 
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respect to the corresponding damage activation function. However, it is also 
necessary to specify how the inactive elastic domain evolves if other damage 
modes are active. It is assumed that the longitudinal and transverse domains 
are not coupled. On the other hand, compression damage is coupled with 
tension damage, as explained in the next section. 


3-4-1 Transverse loading 

As previously described, transverse damage in the form of matrix cracks 
can have different orientations as a result of tension, shear, or compression- 
dominated loads. Under load reversal, transverse cracks, which are perpen- 
dicular to the ply mid-plane, do not affect the compression response: elastic 
domain and the compressive damage variable (d 2 _) are unaffected by r 2+ . 

On the other hand, matrix cracks at a fracture angle cto = 53° caused by high 
compressive transverse stresses have the same effect as cracks perpendicular to 
the mid-plane (a = 0°) when the load is reversed from compression to tension. 
Therefore, the evolution of the transverse tensile clastic domain threshold (r 2+ ) 
is governed by both damage mechanisms. 

Based on the above considerations, the evolution of the elastic domain in the 
transverse direction can be represented by the following equations: 


Tension loading: r 2+ = 0 2+ and r 2 _ = 0 


Compression loading: r 2 _ 


and r 2+ 


<fi 2 - if r 2+ < r 2 _ 
0 if r 2+ > r 2 _ 


The integration of the previous expressions results in: 


r 2+ = max 1 1 
r 2 _ = max { 1 


max 

s=0,t 

max 

s=0,t 


{<«-} 

{•%-} 


. max 

s=0,t 



(24) 


3-4-2 Longitudinal loading 

Under longitudinal tensile stresses, the fracture plane is perpendicular to the 
fiber direction. When reversing the load, the cracks close and can still transfer 
load. However, the broken and misaligned fibers do not carry any additional 


18 



load. Therefore, the compressive stiffness is influenced by longitudinal damage, 
ffowever, the elastic domain is assumed to remain unchanged. 

Under longitudinal compression, damaged material consisting of broken fibers 
and matrix cracks forms a kink band, and there is not a unique orientation 
for the damage planes. When the loads are reversed, the cracks generated in 
compression open and the clastic domain threshold increases. 

Therefore, the evolution of damage thresholds for longitudinal damage are 
defined as: 


Tension loading: ri+ = 4>i+ and r i_ = 0 


Compression loading: fi- 


, . (25) 

! Oi if r L . < r 
and ri + = < 

I 0 if r 1+ > ri- 


The integration of the previous expressions results in: 


r 1+ = max |l, max {^f + } , max } 
ri_ = max 1 1, max 1 1 


(26) 


3.5 Damage evolution laws 


The internal variables tn define the threshold of the elastic domains, and are 
related to the damage state of each lamina, i.e., the damage variables depend 
on the values of the internal variables. In order to fully define the constitutive 
model, it is necessary to define the relation between the internal variables and 
the damage variables. 

When material is undamaged the internal variables rjv take the initial value 
of 1, and djv(rAr = 1) = 0. Equations (24) and (26) define the evolution of the 
internal variables assuring that rjv > 0. As shown in equations (18) and (19), 
the condition for positive dissipation is satisfied if djv > 0. The condition for 
positive dissipation is automatically fulfilled if the damage evolution law satis- 
fies the condition ddpf/dr N > 0. When the material is completely damaged, a 
fracture plane is created, the strains are localized in a plane in which rjv — > oo 
and the related components of the stiffness tensor are zero, djv(r/v — > oo) = 1. 

Matrix cracks are related to the internal variables r 2 + and ri— The internal 
variable r-i- accounts for compressive damage only, whereas ?~ 2 + accounts for 
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both, compressive and tensile, damage. Therefore, for positive transverse nor- 
mal stresses damage is a function of d 2+ (r 2+ ) because both types of cracks 
(a = 0° and a = 53°) are open. Under transverse compressive loads, the 
damage is only influenced by the inclined cracks, c/ 2 _ (r 2 _). 

Kink bands are related to the internal variable ri_. The internal variable r 1+ 
accounts for kink bands and fiber tensile fracture. For positive longitudinal 
normal stresses, the material loses stiffness as a result of both damage modes 
because the cracks open. Therefore, the damage variable can be expressed as 
di+ (n+). 

When a lamina which is fully damaged in tension (di + = 1) is subjected to 
load reversal and the crack closes, some of the original stiffness is recovered 
because the tractions can be transmitted through the crack faces. However, 
the broken fibers lose their alignment. Assuming that the fibers do not carry 
any load, which can be considered as a limit case, the compressive stiffness can 
be approximated by the rule of mixtures applied for components in parallel 
as: (1 — di_)Ei = V m E m . The generalization of the above arguments for an 
intermediate damage state can be expressed as d\_ = A±d\ + , with: 




VfE f 


V m E m + VfE 


f^f 


E\ — E‘2 

) 

E x 


(27) 


where Ef and E m are the fiber and matrix Young modulus, Vf and V m the 
corresponding volume fractions, and b is an adjustment parameter between 0 
and 1. If b — 1 the stiffness recovery is due only by the matrix, and if b = 0, 
the stiffness recovery is total and it is assumed that broken fibers do not lose 
alignment under compressive loads, and that the initial stiffness is recovered. 

Fiber damage (di±) is not influenced by matrix cracking (r 2 ±) as shown in ex- 
perimental results carried by Carlsson and Pipes [50] and in micromechanical 
models [51]- [54] of cracked composites. Therefore, the longitudinal stiffness is 
not function of matrix transverse cracks. 


The shear stiffness is reduced as a result of longitudinal and transverse cracks 
regardless of their orientation. Under these circumstances, the damage variable 
dg is given by: 


d 6 = l-[l-d 6 (r 2+ )](l-d 1+ ) (28) 

The damage evolution laws used force strain-softening as soon as one damage 
activation criterion is satisfied. Softening constitutive equations result in phys- 
ically inadmissable responses: the damage is localized in a plane and fracture 
occurs without energy dissipation. The numerical implementation of soften- 
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ing constitutive equations using the finite element method results in mesh- 
dependent results because the energy dissipated is a function of the element 
size. 

The solution normally used to ensure the correct computation of the energy 
dissipated regardless of the refinement of the mesh is to adjust the damage 
evolution laws using a characteristic dimension of the finite element. The def- 
inition of the damage evolution laws is therefore related to the computational 
model, and will be described in the next section. 


4 COMPUTATIONAL MODEL 


Two different inelastic responses have to be taken into account in the numer- 
ical simulation of damage using Continuum Damage Mechanics. While the 
stress-strain response of a material exhibits a positive-definite tangent stiff- 
ness tensor, the damage zone increases along all the directions. This initial 
stage of damage is commonly referred to as diffuse damage, and the numerical 
solution is independent of the numerical discretization. For example, matrix 
cracking in a multidirectional composite is a form of diffuse damage if the 
kinematics of laminate theory is used. 

When the tangent stiffness tensor is not positive definite, damage localizes 
in a narrow band, and the numerical solution depends upon the numerical 
discretization: decreasing the element size in the localized zone decreases the 
computed energy dissipated. Therefore, the structural response is not objective 
because it does not converge to a unique solution with mesh refinement. 

The proposed damage model uses a constitutive model that forces localiza- 
tion as soon as one of the damage activation functions associated with the 
onset of transverse or longitudinal cracking is satisfied, i.e., when F n = 0. 
In order to guarantee that the numerical solution is independent of the dis- 
cretization a characteristic element length is used in the constitutive model 
using a procedure based on the crack band model proposed by Bazant [26]. 


4-1 Damage laws in softening regime: crack band model 


Bazant’s crack band model [26] assures the objective response of the global 
finite element model by regularizing the computed dissipated energy using a 
characteristic dimension of the finite element and the fracture toughness: 
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, M = l+,l-,2+,2-,6 


(29) 


Gm 


where Gm is the fracture toughness, (Jm is the energy dissipated per unit 
volume, and l* is the characteristic length of the finite element. For square 
elements, with an aspect ratio approximately equal to one, the characteristic 
element length can be approximated by the following expression [26]: 


_ VAip 
cos (7) 

where |y| < 45° is the angle of the mesh lines with the crack direction and Ajp 
is the area associated at each integration point. For an unknown direction of 

— 2L 

crack propagation, the average of this expression can be used, l* — | f 0 4 l*d 7 = 

1 . 12 \/ Ajp. 

When the crack propagation path can be estimated in advance, it is recom- 
mended to align the mesh with the direction of crack propagation because 
cracks tend to evolve along the mesh lines. If the crack propagation is aligned 
with the mesh lines, the characteristic length must be the square root of the 
area corresponding to an element integration point, i.e., 7 = 0. 

For triangular elements, the typical characteristic length is determined by the 
expression: 



r=2 /W (3i) 

A more accurate measure of the characteristic element length would be ob- 
tained using the element projections for both possible crack directions, trans- 
verse and longitudinal. 

The crack band model assumes that the failure process zone can be represented 
by a damaged finite element zone of one element width. This approximated 
method for achieving the objectivity of the global response is appropriate for 
the treatment of large structures under complex damage mechanisms, such as 
the ones occurring in advanced composites. 

The exponential damage evolution laws proposed here are expressed in the 
following general form: 


dM — 1 — 


^77 T exp {Am [1 - f N (rjv)]} / (r K ) 


(32) 
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where the function (rjv) is selected to force the softening of the constitu- 
tive relation. / ( tk ) is the coupling factor between damage laws and elastic 
threshold domains. 

The damage evolution laws for each damage variable are: 


di+ = 1 - ^ exp[A 1+ (l - r 1+ )] 

<l\- - I - n 'd exp[Ai_(l - n-)]f(Af, r i + ) 

d'2+ = 1 - h+ lr 2+ ) exp[A 2+ (l - / 2 +(r 2+ ))] (33) 

d 2 ~ = 1 - exp[A 2 _(l - r 2 _)] 

d 6 (r 2+ ) = 1 - ^ exp[A 6 (l - r 2+ )] 

where / 2+ (r 2+ ) is a function of the same order as the damage onset function 
in order to force softening: 


/ 2+ (r 2 +) = Y g { 9 ~ 1 + “ df + 4 3 r 2+) (34) 

The coupling factor Af is used for the interaction of elastic domains in the 
longitudinal (fiber) direction, and it is defined from equation d \ _ = Afd\ + 
and (33) for ri_ = 1 as: 


/ (4dn) = 1 - Af + Af — exp [A 1+ (1 - ri+)] (35) 

r i+ 

where Af is a material parameter defined in equation (27). 

The energy dissipated per unit volume for uniaxial stress conditions is obtained 
by integrating the rate of dissipation, equation (18): 

ro o . roc QQ Od'A/r 

9m — / Y M dMdt= / — — — — -dr M: M— 1+, 1— , 2+, 2— , 6 (36) 

Jo J i od M or M 

Applying the crack band model, equation (29): 


f°° d G dd M , Gm 
J i dd M dr M l* 


M 


1 +, 1 - 2 +, 2-,6 


(37) 
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Using equations (29) and (36) and substituting in (37), it is possible to cal- 
culate the damage law parameters Am that assure that the dissipated energy 
computed by the numerical model is independent of mesh refinement, ft is 
possible to obtain analytical closed form solutions for two of the equations 
(37): 


A 1+ — 


21 * 

2E l G l+ - l* X\ 


(38) 


21* S 2 L 

2 G 12 G 6 - l* SI 


(39) 


The remaining parameters, A 2 ± and Ai_, are calculated numerically using 
the algorithm presented in Appendixes A and B. The adjusting parameters 
complete the definition of the constitutive model. 

The material response under load reversal cycles in the transverse and longitu- 
dinal direction resulting from the proposed constitutive model are illustrated 
in Figures 9 and 10 respectively. 



Fig. 9. Transverse load cycle: O-A-B-O-C-D-O-E. 
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Fig. 10. Longitudinal load cycle: O-A-B-O-C-D-E-O-F. 

In the transverse load cycle shown in Figure 9, it can be seen that damage 
created under tensile stresses (O-A-B-O) does not influence the compressive 
behavior (O-C), the size of the elastic domain (r 2 _ = 1), or the damage 
variable (d 2 - = 0). On the other hand, damage created under compressive 
loading (C-D-O) increases the elastic domain in tension (O-E) and affects the 
damage variable d 2+ . 

In the longitudinal load cycle shown in Figure 10, it can be seen that damage 
created under tensile loads (O-A-B-O) influences the compressive behavior 
under load reversal (O-C), corresponding to a decrease of stiffness due the 
misalignment of the fibers (di_ = A^d\ + ). However, the clastic domain size 
remains unchanged (r'i_ = 1). 

For damage created under compressive stresses, two regions can be distin- 
guished. For rq_ < ri + (C-D), the tensile elastic domain threshold does not 
change. For ri_ > r 1+ (D-E), both elastic domains thresholds increase to 
reflect the fracture of fibers in compression. 


4-1.1 Critical finite element size 

The constitutive model must not lead to a local snap-back in the stress-strain 
relation. In other words, the elastic energy of an element at the onset of local- 
ization, which is X\ 4 (l*) 2 1/ (2 Em) with M = 1±,2±,6, must be lower than 
or equal to the fracture energy, GA/Cf. 
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Therefore, the maximum size for the finite element for each damage law M is: 


r < 2E ^ M ,M = 1±, 2±, 6 (40) 

X M 

where Em, Gm and Xm are the Young modulus, fracture energies and strengths, 
respectively. 

If a finite element model consist of elements larger that the maximum size 
prescribed by equation (40), there is an alternative to further refining the 
mesh. The snap-back in the constitutive model can be avoided by reducing 
the corresponding strength [26] while taking the parameter Am to infinity: 


Am = (41) 

Under these circumstances, the damage variables take two possible values: 
d,M = 0 if tty = 1 or oUf = 1 if rjv > 1. It is clear that the modification of the 
strengths, although assuring the correct calculation of the energy dissipated, 
should not be performed in the elements representing the region where crack 
initiation, which is controlled by the stress tensor, takes place. In the regions 
of stress concentrations, where crack initiation is likely to take place, the 
mesh should be sufficiently refined to avoid any adjustment of the material 
properties. 


4-1.2 Fracture toughness 

Each damage evolution function includes one adjusting parameter, Am, M = 
1±, 2±, 6, that needs to be calculated using the corresponding component of 
the fracture toughness, Gm, representing the energy dissipated by inelastic 
processes in the fracture process zone. 

G- 2 + and G 6 correspond to the fracture toughness of a transverse crack in 
mode I and II, respectively. The mode I component of the fracture tough- 
ness, C? 2 +, can be measured using the Double Cantilever Beam (DCB) test 
(ASTM-D5528). The mode II component of the fracture toughness, Gq, can 
be measured using the Four-Point End Notched Flexure (4-ENF) test speci- 
men [28]. 

G l+ corresponds to the mode I component of the fracture toughness for a 
longitudinal crack. There is no standard test method to measure this property. 
The suggested test method to measure G\+ is the Compact Tension (CT) test 
specimen proposed by Pinho et al. [28] 
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The measurement of the energy dissipated that is associated with longitudinal 
compressive loading is far more complex because several complex dissipative 
phenomena are involved, such as crack growth, crushing and friction. These 
events occur sequentially throughout the fracture process and transverse in- 
clined cracks and longitudinal kink bands are the result of the combination of 
these dissipative mechanisms. 

Bazant et al. [55] proposed the following expression to evaluate the energy 
dissipated in a kink band, G\-\ 


G\- = W G 6 (42) 

s 

where w is the kink band thickness and s is the distance between two ma- 
trix cracks. This approximation requires a good knowledge of the kink band 
geometry in advance, which is a function of the external loads, the geometry 
of structure and the thickness confinement [56]. This approximation does not 
take into account other dissipative mechanisms in the material such as slid- 
ing of the crack faces. An alternative procedure to measure G i_ is to use the 
compact compression (CC) specimen that triggers kink bands in laminated 
composites, as proposed by Pinho [28]. 

The fracture toughness for transverse compression loading, G 2 - can be calcu- 
lated approximately using the mode II component of the fracture toughness, 
the fracture angle ao and a term accounting for friction between the crack 
faces: 


G‘ 2 - = 6 b at/j,Yc cos ao ~ - b atrj T Yc cos ao (43) 

cos ao cos ao 

where ao ~ 53 ± 3°, t is the lamina thickness, and a is an adjustment para- 
meter between 0 (in an unidirectional laminate) and 1 (in a strongly confined 
lamina) . 


4-1.3 Viscous regularization 

Strain-softening constitutive models cause convergence difficulties when using 
global solution methods, especially for damage in the longitudinal (fiber) di- 
rection. In order to improve the convergence of the numerical algorithm, an 
artificial Duvaut-Lions viscosity model [57] is implemented. The time deriv- 
atives of the internal variables associated with longitudinal failure can be 
defined as: 
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V 


(44) 


r 1+ = 


(max 


S 

l+» ' 


<l>i-}-ri+) 


and ri_ = 


*!- 


where r/ is the viscous parameter. When rj tends to zero and (max |0f + , j — 
r 1+ ) > 0, the mathematical definition of a derivative is obtained and the 
functions (44) tend to the damage thresholds evolution functions defined by 
equations (25). 


A numerical algorithm needs to be implemented for the time integration of 
the internal variables. Using a backward- Euler scheme, the internal variables 
can be updated as: 


r n+l 
' l- 


r n + 1 
r l+ 


max 

max 


i-) 


V 

T) + At 

,n + 1 


r n_ + 


At 

77 + At 


b n+l 


... V r n | At 

1+ ’ 7 ri + At 1+ r) + At 



(45) 


Although some materials exhibit time-dependent response, this regularization 
is implemented with the objective of improving the numerical convergence of 
the model. 

Two undesirable consequences can occur when increasing the viscous para- 
meter. The first one is that the tangent relation becomes positive definite at 
first stage of damage therefore localization is not ensured at damage onset. 
Secondly, the energy dissipated at a material integration point undergoing 
damage evolution increases with the viscous parameter. 


4-2 Material tangent constitutive tensor and algorithm 


The fast convergence rate of the solution algorithm for the non-linear problem 
requires the correct computation of the material tangent constitutive tensor, 

C T : 


a = Ct : £ 


(46) 


where: 


C T = H' 1 : (I - M) 


(47) 



H 1 is the secant constitutive tensor, I is the identity tensor and the tensor 
M is defined as: 


M 


£11 


ddi 


dd\ 


ddi 


(l-d 1 j 2 E 1 den 

{1-dxfEx d£22 

(1-dxfEx de 12 

(722 dd2 

(722 dd2 

(722 dd2 

(1 -d 2 ) 2 E 2 (ten 

(1 -d 2 ) 2 E 2 de 22 

(1 -d 2 ) 2 E 2 dei 2 

<712 dd R 

t -\ .1 \2 /~i . - 

<7i2 ddg 

/-i .i \2 £1. 

<7 12 dda 

/- 1 .7 \2/~l £1. 


(48) 


The scalar components of the tensor M are presented in Appendix C. 

The integration of the constitutive model is performed according to the fol- 
lowing algorithm: 


1 - 

Read the strain tensor at time t 

e* 

2 - 

Compute the effective stress tensor 

a 1 = Hq 1 : e t 

3 - 

Compute the loading functions 

^m(^) 

4 - 

Compute the threshold values 

r Vf( r M > 0 m ) 

5 - 

Compute the damage variables 


6 - 

Compute the nominal stress tensor 

<7* = (H*) -1 : e t 

7 - 

Compute the tangent constitutive tensor 

Q 

II 

1 


It should be noted that the numerical calculation of the adjusting parameters 
Am is only performed once, and only if the damage variable at the integration 
point is greater than zero. Also, the verification of the condition Am > 0 
is performed in the beginning of the analysis, and, if it is not satisfied, the 
material strengths are reduced according to the equation (41). 


5 VALIDATION 


The model developed was implemented in ABAQUS non-linear finite element 
code [58] using a user-subroutine UMAT. 

The model is validated by comparing the predicted failure loads of quasi- 
isotropic laminates containing a central hole and loaded in tension with the 
corresponding experimental data. 
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The experimental data used was previously obtained by Tan [59] for the 
open-hole tension test specimen shown in Figure 11. The material consists 
of T300/1034-C carbon fiber reinforced epoxy with a nominal ply thickness of 
0.1308 mm. The elastic properties and unidirectional strengths are shown in 
Tables 1 and 2 respectively. 
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Fig. 11. Configuration of the open-hole tension test specimen. 


Table 1 

T300/1034-C elastic properties [35], [59]. 


E x (GPa) 

E 2 (GPa) 

G \ 2 (GPa) u \ 2 

146.8 

11.4 

6.1 0.30 


Table 2 

T300/1034-C unidirectional strengths (MPa) [59], [35]. 


Xt 


Yt 

Yc 

Sl 

1730.0 

1379.0 

66.5 

268.2 

58.7 


The components of the fracture toughness for the different damage models 
are required for both the determination of the in-situ strengths and of the 
adjusting parameters A^. 
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The components of the fracture toughness associated with matrix cracking 
used in the model were measured by Shahid and Chang [60] for T300/934 
carbon-epoxy laminates. The components of the fracture toughness associ- 
ated with longitudinal tensile and compressive fracture used were obtained by 
Pinho [28] for a carbon-epoxy composite using the same fiber type (T300/913). 
The values of fracture toughness used in the model are shown in Table 3. 

Table 3 

Fracture toughness (N/mni). 


g 2+ 

g 6 

G 2 - 

Gi+ 

G\- 

0.23 

0.46 

0.76 

89.83 

78.27 


The damage activation functions are established in terms of the in-situ strengths 
The in-situ strengths are calculated from the closed-form equations proposed 
in [36], using a shear response factor [3 = 3.2 x ICC 8 mm 6 /N 3 . The resulting 
in-situ strengths are shown in Table 4. 

Table 4 

In-situ strengths (MPa). 



Embedded ply 

Outer ply 

Thick ply 

Y t 

158.8 

101.2 

105.4 

Sl 

109.5 

89.8 

73.4 


The coefficients of thermal expansion of the material used are an = —1.0 x 
10 _6 /°C and a 2 2 = 26.0 x 10“ 6 /°C, and the temperature change is AT = 
— 152°C [61]. The effects of moisture absorption are neglected in the model. 

The lay-ups tested in specimens with the geometry shown in Figure 11 are 
[07[±45 o ] 3 /90^] s ,[07[±45 o ] 2 /90^ and [07 ± 45790?],. 

The finite element model created use 4-node shell elements. Based on the prop- 
erties reported in Tables 1-4, the maximum element size that avoids lowering 
the strength is 0.508 mm. The mesh in the vicinity of the hole, corresponding 
to the region where damage occurs, consist of elements with sizes ranging from 
0.127 mm. (hole edge) to 0.635 mm. (specimen edge). 

It is worth noting that shell elements based on lamination theory assume 
a linear strain field along laminate thickness expressed by two tensors: mid 
plane strains (e 0 ) and curvatures (k). This simplified kinematic description 
may not be able to detect the localization of strains in a ply constrained by 
sub-laminates. In other words, damage mechanisms, such as matrix trans- 
verse cracking in just one ply of a multidirectional laminate, correspond to 
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distributed or diffuse damage if lamination theory is employed. Localization 
of damage occurs when the different damage mechanisms occur in all plies of 
the laminate. 

Before the localization of damage in the laminate, the ply constitutive models 
that can be used are based on strain softening, such as in the model proposed 
here, or on elastic analysis of cracked plies. The latter solution is normally only 
obtainable for a periodic distribution of transverse matrix cracks in central 90° 
plies of rectangular laminates under constant stresses [53], [54], 

Figure 12 shows the load-displacement relation of the three specimens simu- 
lated. 


TJ 

cS 

O 

hJ 



0,0 0,2 0,4 0,6 0,8 1,0 1,2 

Applied end displacement (mm) 


Fig. 12. Predicted load-displacement relations. 

Table 5 presents the predicted and experimental failure stresses, a N , defined 
using the failure load P u , and the specimen width and thickness, w and t 
respectively, as a N = ^ = 9.71 P u . 

Table 5 

Predicted and measured failure stress, a N (MPa). 


Lay-up 

Experimental [59] 

Predicted 

Error (%) 

[0°/[±45 o ] 3 /90^ 

235.8 

225.5 

-4.4 

[0°/[±45°] 2 /90^] s 

185.8 

192.4 

3.7 

[0°/ ± 45°/90?] s 

160.0 

139.3 

-12.9 


The predicted failure loads are in excellent agreement with the experimental 
failure loads measured by Tan [59]. 
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6 CONCLUSIONS 


A new constitutive model for the prediction of damage onset, growth and 
ultimate failure of composite structure under plane stress was proposed. The 
onset of the different intralaminar damage mechanisms is predicted using a 
simplification of the LaRC04 failure criteria. 

The constitutive model proposed is based on four ply fracture planes and 
accounts for the unilaterality of damage by its ability to represent complex 
load histories, including tension-compression load reversals. 

The constitutive law was implemented in a computational model that en- 
sures that the computed dissipated energy is independent of the discretization. 
Therefore, the numerical solution is objective with respect to mesh refinement. 

The computational model developed was used in the simulation of open-hole 
test specimens loaded in tension using different lay-ups. An excellent agree- 
ment between the predicted and measured failure stresses was obtained. 
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Appendix A: Integration Algorithm 


In order to calculate the adjusting parameters Am used in the damage evolu- 
tion laws, it is necessary to integrate the following equation numerically: 


9m 



1 -^ 12^21 V d 

(1 — (Im) ^ 12 ^ 21 ) 2 E m dr N 


(A-l) 


The Simpson method of numerical integration approximates the solution using 
quadratic polynomials. The general form of the polynomial can be expressed 
as: 
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where h is the step increment, and f l M = 
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Since the damage laws selected tend to zero, it is necessary to define a point to 
stop the integration. When the stress becomes less than K times the stress at 
the onset of localization, the remaining energy can be neglected. The increment 
(h) can be selected defining the number of steps (n) as: 


hi 



(A-3) 


The algorithm is implemented as: 
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1- Select the parameters n and K 


2- Initialize 

3- Compute step size 

4- WHILE CONT<n 


r=l , G=0 and CONT=0 
h 


DO 1=1:3 


f(I) = 


I — U12U21 

1 — (1— d M )t/12^21 


°"m dd M 
2 E m dr N 


r=r+h 


END DO 


r=r-h 

9 = 9 + | (/ (1) + 4/ (2) + / (3)) 
CONT=CONT+l 


END WHILE 




Appendix B: Secant method to determine the parameters Am 


To find the adjustment parameters for the damage law, it is necessary to 
integrate the stress-strain relation in terms of the unknown parameter Am- The 
integration can be done numerically with the algorithm presented in Appendix 
A. To iterate on the value of A M , the secant method is used. The problem to 
be solved can be expressed as: 

9m (Am) = 0 (B-l) 

To select the two parameters to start the iteration, the following approxima- 
tion is used: 


07* v2 
^ -A M 


/* X 2 

1 ^ M 


and A° m = 0.5 A l M 


The function g l (Aj) is only defined for positives values of A. Defining the 
minimization function as: 



(B-3) 


the following algorithm is proposed: 
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Appendix C: Material tangent stiffness tensor 


If artificial viscosity is considered, incremental damage laws have to be imple- 
mented: 
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Applying the chain rule, the damage evolution can be written as: 
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Derivation of softening damage laws: 
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Evolution of damage thresholds respect the strains: 
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